home *** CD-ROM | disk | FTP | other *** search
/ Night Owl 6 / Night Owl's Shareware - PDSI-006 - Night Owl Corp (1990).iso / 039a / jpsrc2.zip / JFWDDCT.C < prev    next >
C/C++ Source or Header  |  1991-12-02  |  7KB  |  227 lines
  1. /*
  2.  * jfwddct.c
  3.  *
  4.  * Copyright (C) 1991, Thomas G. Lane.
  5.  * This file is part of the Independent JPEG Group's software.
  6.  * For conditions of distribution and use, see the accompanying README file.
  7.  *
  8.  * This file contains the basic DCT (Discrete Cosine Transform)
  9.  * transformation subroutine.
  10.  *
  11.  * This implementation is based on Appendix A.2 of the book
  12.  * "Discrete Cosine Transform---Algorithms, Advantages, Applications"
  13.  * by K.R. Rao and P. Yip  (Academic Press, Inc, London, 1990).
  14.  * It uses scaled fixed-point arithmetic instead of floating point.
  15.  */
  16.  
  17. #include "jinclude.h"
  18.  
  19.  
  20. /* We assume that right shift corresponds to signed division by 2 with
  21.  * rounding towards minus infinity.  This is correct for typical "arithmetic
  22.  * shift" instructions that shift in copies of the sign bit.  But some
  23.  * C compilers implement >> with an unsigned shift.  For these machines you
  24.  * must define RIGHT_SHIFT_IS_UNSIGNED.
  25.  * RIGHT_SHIFT provides a signed right shift of an INT32 quantity.
  26.  * It is only applied with constant shift counts.
  27.  */
  28.  
  29. #ifdef RIGHT_SHIFT_IS_UNSIGNED
  30. #define SHIFT_TEMPS    INT32 shift_temp;
  31. #define RIGHT_SHIFT(x,shft)  \
  32.     ((shift_temp = (x)) < 0 ? \
  33.      (shift_temp >> (shft)) | ((~0) << (32-(shft))) : \
  34.      (shift_temp >> (shft)))
  35. #else
  36. #define SHIFT_TEMPS
  37. #define RIGHT_SHIFT(x,shft)    ((x) >> (shft))
  38. #endif
  39.  
  40.  
  41. /* The poop on this scaling stuff is as follows:
  42.  *
  43.  * We have to do addition and subtraction of the integer inputs, which
  44.  * is no problem, and multiplication by fractional constants, which is
  45.  * a problem to do in integer arithmetic.  We multiply all the constants
  46.  * by DCT_SCALE and convert them to integer constants (thus retaining
  47.  * LG2_DCT_SCALE bits of precision in the constants).  After doing a
  48.  * multiplication we have to divide the product by DCT_SCALE, with proper
  49.  * rounding, to produce the correct output.  The division can be implemented
  50.  * cheaply as a right shift of LG2_DCT_SCALE bits.  The DCT equations also
  51.  * specify an additional division by 2 on the final outputs; this can be
  52.  * folded into the right-shift by shifting one more bit (see UNFIXH).
  53.  *
  54.  * If you are planning to recode this in assembler, you might want to set
  55.  * LG2_DCT_SCALE to 15.  This loses a bit of precision, but then all the
  56.  * multiplications are between 16-bit quantities (given 8-bit JSAMPLEs!)
  57.  * so you could use a signed 16x16=>32 bit multiply instruction instead of
  58.  * full 32x32 multiply.  Unfortunately there's no way to describe such a
  59.  * multiply portably in C, so we've gone for the extra bit of accuracy here.
  60.  */
  61.  
  62. #ifdef EIGHT_BIT_SAMPLES
  63. #define LG2_DCT_SCALE 16
  64. #else
  65. #define LG2_DCT_SCALE 15    /* lose a little precision to avoid overflow */
  66. #endif
  67.  
  68. #define ONE    ((INT32) 1)
  69.  
  70. #define DCT_SCALE (ONE << LG2_DCT_SCALE)
  71.  
  72. /* In some places we shift the inputs left by a couple more bits, */
  73. /* so that they can be added to fractional results without too much */
  74. /* loss of precision. */
  75. #define LG2_OVERSCALE 2
  76. #define OVERSCALE  (ONE << LG2_OVERSCALE)
  77. #define OVERSHIFT(x)  ((x) <<= LG2_OVERSCALE)
  78.  
  79. /* Scale a fractional constant by DCT_SCALE */
  80. #define FIX(x)    ((INT32) ((x) * DCT_SCALE + 0.5))
  81.  
  82. /* Scale a fractional constant by DCT_SCALE/OVERSCALE */
  83. /* Such a constant can be multiplied with an overscaled input */
  84. /* to produce something that's scaled by DCT_SCALE */
  85. #define FIXO(x)  ((INT32) ((x) * DCT_SCALE / OVERSCALE + 0.5))
  86.  
  87. /* Descale and correctly round a value that's scaled by DCT_SCALE */
  88. #define UNFIX(x)   RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1)), LG2_DCT_SCALE)
  89.  
  90. /* Same with an additional division by 2, ie, correctly rounded UNFIX(x/2) */
  91. #define UNFIXH(x)  RIGHT_SHIFT((x) + (ONE << LG2_DCT_SCALE), LG2_DCT_SCALE+1)
  92.  
  93. /* Take a value scaled by DCT_SCALE and round to integer scaled by OVERSCALE */
  94. #define UNFIXO(x)  RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1-LG2_OVERSCALE)),\
  95.                    LG2_DCT_SCALE-LG2_OVERSCALE)
  96.  
  97. /* Here are the constants we need */
  98. /* SIN_i_j is sine of i*pi/j, scaled by DCT_SCALE */
  99. /* COS_i_j is cosine of i*pi/j, scaled by DCT_SCALE */
  100.  
  101. #define SIN_1_4 FIX(0.707106781)
  102. #define COS_1_4 SIN_1_4
  103.  
  104. #define SIN_1_8 FIX(0.382683432)
  105. #define COS_1_8 FIX(0.923879533)
  106. #define SIN_3_8 COS_1_8
  107. #define COS_3_8 SIN_1_8
  108.  
  109. #define SIN_1_16 FIX(0.195090322)
  110. #define COS_1_16 FIX(0.980785280)
  111. #define SIN_7_16 COS_1_16
  112. #define COS_7_16 SIN_1_16
  113.  
  114. #define SIN_3_16 FIX(0.555570233)
  115. #define COS_3_16 FIX(0.831469612)
  116. #define SIN_5_16 COS_3_16
  117. #define COS_5_16 SIN_3_16
  118.  
  119. /* OSIN_i_j is sine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
  120. /* OCOS_i_j is cosine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
  121.  
  122. #define OSIN_1_4 FIXO(0.707106781)
  123. #define OCOS_1_4 OSIN_1_4
  124.  
  125. #define OSIN_1_8 FIXO(0.382683432)
  126. #define OCOS_1_8 FIXO(0.923879533)
  127. #define OSIN_3_8 OCOS_1_8
  128. #define OCOS_3_8 OSIN_1_8
  129.  
  130. #define OSIN_1_16 FIXO(0.195090322)
  131. #define OCOS_1_16 FIXO(0.980785280)
  132. #define OSIN_7_16 OCOS_1_16
  133. #define OCOS_7_16 OSIN_1_16
  134.  
  135. #define OSIN_3_16 FIXO(0.555570233)
  136. #define OCOS_3_16 FIXO(0.831469612)
  137. #define OSIN_5_16 OCOS_3_16
  138. #define OCOS_5_16 OSIN_3_16
  139.  
  140.  
  141. /*
  142.  * Perform a 1-dimensional DCT.
  143.  * Note that this code is specialized to the case DCTSIZE = 8.
  144.  */
  145.  
  146. INLINE
  147. LOCAL void
  148. fast_dct_8 (DCTELEM *in, int stride)
  149. {
  150.   /* many tmps have nonoverlapping lifetime -- flashy register colourers
  151.    * should be able to do this lot very well
  152.    */
  153.   INT32 in0, in1, in2, in3, in4, in5, in6, in7;
  154.   INT32 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
  155.   INT32 tmp10, tmp11, tmp12, tmp13;
  156.   INT32 tmp14, tmp15, tmp16, tmp17;
  157.   INT32 tmp25, tmp26;
  158.   SHIFT_TEMPS
  159.   
  160.   in0 = in[       0];
  161.   in1 = in[stride  ];
  162.   in2 = in[stride*2];
  163.   in3 = in[stride*3];
  164.   in4 = in[stride*4];
  165.   in5 = in[stride*5];
  166.   in6 = in[stride*6];
  167.   in7 = in[stride*7];
  168.   
  169.   tmp0 = in7 + in0;
  170.   tmp1 = in6 + in1;
  171.   tmp2 = in5 + in2;
  172.   tmp3 = in4 + in3;
  173.   tmp4 = in3 - in4;
  174.   tmp5 = in2 - in5;
  175.   tmp6 = in1 - in6;
  176.   tmp7 = in0 - in7;
  177.   
  178.   tmp10 = tmp3 + tmp0;
  179.   tmp11 = tmp2 + tmp1;
  180.   tmp12 = tmp1 - tmp2;
  181.   tmp13 = tmp0 - tmp3;
  182.   
  183.   in[       0] = (DCTELEM) UNFIXH((tmp10 + tmp11) * SIN_1_4);
  184.   in[stride*4] = (DCTELEM) UNFIXH((tmp10 - tmp11) * COS_1_4);
  185.   
  186.   in[stride*2] = (DCTELEM) UNFIXH(tmp13*COS_1_8 + tmp12*SIN_1_8);
  187.   in[stride*6] = (DCTELEM) UNFIXH(tmp13*SIN_1_8 - tmp12*COS_1_8);
  188.  
  189.   tmp16 = UNFIXO((tmp6 + tmp5) * SIN_1_4);
  190.   tmp15 = UNFIXO((tmp6 - tmp5) * COS_1_4);
  191.  
  192.   OVERSHIFT(tmp4);
  193.   OVERSHIFT(tmp7);
  194.  
  195.   /* tmp4, tmp7, tmp15, tmp16 are overscaled by OVERSCALE */
  196.  
  197.   tmp14 = tmp4 + tmp15;
  198.   tmp25 = tmp4 - tmp15;
  199.   tmp26 = tmp7 - tmp16;
  200.   tmp17 = tmp7 + tmp16;
  201.   
  202.   in[stride  ] = (DCTELEM) UNFIXH(tmp17*OCOS_1_16 + tmp14*OSIN_1_16);
  203.   in[stride*7] = (DCTELEM) UNFIXH(tmp17*OCOS_7_16 - tmp14*OSIN_7_16);
  204.   in[stride*5] = (DCTELEM) UNFIXH(tmp26*OCOS_5_16 + tmp25*OSIN_5_16);
  205.   in[stride*3] = (DCTELEM) UNFIXH(tmp26*OCOS_3_16 - tmp25*OSIN_3_16);
  206. }
  207.  
  208.  
  209. /*
  210.  * Perform the forward DCT on one block of samples.
  211.  *
  212.  * A 2-D DCT can be done by 1-D DCT on each row
  213.  * followed by 1-D DCT on each column.
  214.  */
  215.  
  216. GLOBAL void
  217. j_fwd_dct (DCTBLOCK data)
  218. {
  219.   int i;
  220.   
  221.   for (i = 0; i < DCTSIZE; i++)
  222.     fast_dct_8(data+i*DCTSIZE, 1);
  223.  
  224.   for (i = 0; i < DCTSIZE; i++)
  225.     fast_dct_8(data+i, DCTSIZE);
  226. }
  227.